/**
   \file  image_utils.hh
   \brief Image processing utility routines.

   This file provides some utility functions that work with Boost.GIL
   images.
*/

/*
   This file is part of libgist.

   libgist is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2 of the License, or (at your
   option) any later version.

   libgist is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public License
   along with libgist; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
*/

/*
   REVISION HISTORY

   $HeadURL: http://libgist.googlecode.com/svn/branches/siagian_dev/src/image_utils.hh $
   $Id: image_utils.hh 60 2011-09-11 07:33:47Z libgist@gmail.com $
*/

#ifndef LIBGIST_IMAGE_UTILS_DOT_HH
#define LIBGIST_IMAGE_UTILS_DOT_HH

//------------------------------ HEADERS --------------------------------

// Boost.GIL
#include <boost/gil/color_base_algorithm.hpp>
#include <boost/gil/image.hpp>
#include <boost/gil/utilities.hpp>

// Boost.Lambda
#include <boost/lambda/lambda.hpp>

// Standard C++
#include <algorithm>
#include <vector>

// Standard C
#include <math.h>

//----------------------- NAMESPACE DEFINITION --------------------------

namespace gist {

//--------------------------- INTERPOLATION -----------------------------

/**
   \brief Quadruple input image size, filling in extra pixels using
   linear interpolation (i.e., averaging) of neighbouring pixels from
   source image.

   \param  input_image The input image that is to be interpolated.
   \return The interpolated image.

   The interpolation process may be described with the following diagram
   showing the layout of pixels in the output image:

   \verbatim
                                 a  A  b
                                 B  C  D
                                 c  E  d
   \endverbatim

   The above diagram shows a small patch of the newly created image. The
   lower-case letters indicate pixels from the original image whereas the
   upper-case letters are interpolated pixels. The interpolated pixels
   are computed as shown below:

   \verbatim
                         A = (a + b) / 2
                         B = (a + c) / 2
                         C = (a + b + c + d) / 4
                         D = (b + d) / 2
                         E = (c + d) / 2
   \endverbatim

   The last row and column are processed specially using duplication.

   NOTE: As this is a template function, it can work with any Boost.GIL
   image types. Note also that we specify the output image type first and
   then the input image type in the template's type list. This allows
   clients to get away with specifying only the output image type and
   relying on type deduction for the input image type.
*/
template<typename output_image_type, typename input_image_type>
output_image_type interpolate(const input_image_type& input_image)
{
   namespace gil = boost::gil ;

   typedef typename input_image_type::const_view_t in_view ;
   typedef typename in_view::xy_locator in_loc ;
   typedef typename in_loc::cached_location_t in_loc_cached ;
   typedef typename output_image_type::view_t out_view ;
   typedef typename out_view::xy_locator out_loc ;
   typedef typename out_loc::cached_location_t out_loc_cached ;

   in_view in = gil::const_view(input_image) ;
   in_loc src = in.xy_at(0,0) ; // source pixel
   in_loc_cached src_e  = src.cache_location(1,0) ; // right of source pixel
   in_loc_cached src_s  = src.cache_location(0,1) ; // below source pixel
   in_loc_cached src_se = src.cache_location(1,1) ; // diagonally right below

   output_image_type output_image(in.width() * 2, in.height() * 2) ;
   out_view out = gil::view(output_image) ;
   out_loc  dst = out.xy_at(0,0) ; // dest. loc. corresponding to source pixel
   out_loc_cached dst_e  = dst.cache_location(1,0) ; // right of dest. pixel
   out_loc_cached dst_s  = dst.cache_location(0,1) ; // below dest. pixel
   out_loc_cached dst_se = dst.cache_location(1,1) ; // diagonally right below

   using boost::lambda::_1 ;
   using boost::lambda::_2 ;
   using gil::static_transform ;

   // Interpolate all rows except the last one
   for (int y = 0; y < in.height() - 1; ++y)
   {
      // Interpolate all columns on current row except last one
      for (int x = 0; x < in.width() - 1; ++x, ++src.x(), ++dst.x(), ++dst.x())
      {
         *dst = *src ; // pixel a (see comment preceding function body)

         // Pixels A and B (see comment preceding function body)
         static_transform(*src, src[src_e], dst[dst_e], (_1 + _2)/2) ;
         static_transform(*src, src[src_s], dst[dst_s], (_1 + _2)/2) ;

         // Pixel C (see comment preceding function body)
         static_transform(*src, src[src_e], dst[dst_se],(_1 + _2)/4) ;
         static_transform(src[src_s],  dst[dst_se], dst[dst_se], _1/4 + _2) ;
         static_transform(src[src_se], dst[dst_se], dst[dst_se], _1/4 + _2) ;
      }

      // Duplicate last column of current row
      *dst = *src ; // pixel a
      dst[dst_e]  = *src ; // pixel A (border)
      static_transform(*src, src[src_s], dst[dst_s], (_1 + _2)/2) ; // pixel B
      dst[dst_se] = dst[dst_s] ; // pixel C (border)

      // Setup for next row
      src += gil::point2<std::ptrdiff_t>(- in.width() + 1, 1) ;
      dst += gil::point2<std::ptrdiff_t>(-out.width() + 2, 2) ;
   }

   // Duplicate last row except last pixel
   for (int x = 0; x < in.width() - 1; ++x, ++src.x(), ++dst.x(), ++dst.x())
   {
      *dst = *src ; // pixel a
      static_transform(*src, src[src_e], dst[dst_e], (_1 + _2)/2) ; // pixel A
      dst[dst_s]  = *src ; // pixel B (border)
      dst[dst_se] = dst[dst_e] ; // pixel C (border)
   }

   // Finally, duplicate bottom-right corner
   *dst = *src ; // pixel a
   dst[dst_e]  = *src ; // pixel A (right border)
   dst[dst_s]  = *src ; // pixel B (bottom border)
   dst[dst_se] = *src ; // pixel C (bottom right corner)

   return output_image ;
}

/**
   \brief Interpolation function when output and input image types are
   the same.

   \param  input_image The input image that is to be interpolated.
   \return The interpolated image.

    Often, the input and output image types are the same. This is a
    convenience function that allows clients to call the generic
    interpolate function with minimum fuss (using parameter type
    deduction) so that the call site does not have to be cluttered with
    unnecessary template verbosity.
*/
template<typename image_type>
inline image_type interpolate(const image_type& input_image)
{
   return interpolate<image_type, image_type>(input_image) ;
}

//---------------------------- CONVOLUTION ------------------------------

/**
   \brief  Convolve input image with supplied kernel along x-direction.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   This function applies the one-dimensional kernel (specified via an
   STL vector of floats) to the input image along the x-direction and
   returns a new image object containing the filteration result.

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename output_image_type, typename input_image_type>
output_image_type
xfilter(const input_image_type& img, const std::vector<float>& kernel)
{
   namespace gil = boost::gil ;

   const int W  = img.width() ;
   const int H  = img.height();
   const int w  = kernel.size() ;
   const int w2 = (w - 1)/2 ;

   // Flip the kernel (this makes the convolution a little more efficient)
   // DEVNOTE: lenrek = kernel spelled backwards, i.e., flipped kernel.
   std::vector<float> lenrek(w) ;
   float kernel_sum = 0 ;
   for (int i = 0; i < w; ++i) {
      kernel_sum  += kernel[i] ;
      lenrek[w - 1 - i] = kernel[i] ;
   }

   typedef typename input_image_type::const_view_t in_view ;
   typedef typename in_view::x_iterator in_iterator ;
   typedef typename in_view::value_type pixel ;

   typedef typename output_image_type::view_t out_view ;
   typedef typename out_view::x_iterator out_iterator ;

   // Create a result image and setup "pointers" into the source and
   // destination pixel arrays.
   output_image_type output_img(W, H) ;
   in_view  in  = gil::const_view(img) ;
   out_view out = gil::view(output_img);

   using gil::static_transform ;
   using boost::lambda::_1 ;
   using boost::lambda::_2 ;

   // Do the convolution...
   if (W < w) // special case for small images (i.e., smaller than the kernel)
   {
      for (int y = 0; y < H; ++y)
      {
         in_iterator  src =  in.row_begin(y) ;
         out_iterator dst = out.row_begin(y) ;
         for (int x = 0; x < W; ++x)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = 0; k < w; ++k)
               if (x + k - w2 >= 0 && x + k - w2 < W) {
                  static_transform(c, src[x + k - w2], c, _1 + _2*lenrek[k]) ;
                  n += lenrek[k] ;
               }
            static_transform(c, dst[x], _1 * kernel_sum/n) ;
         }
      }
   }
   else // convolution of regular sized images (i.e., larger than the kernel)
   {
      for (int y = 0; y < H; ++y)
      {
         in_iterator  src =  in.row_begin(y) ;
         out_iterator dst = out.row_begin(y) ;

         // Left-most w2 pixels
         for (int x = w2; x > 0; --x, ++dst)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = x; k < w; ++k) {
               static_transform(c, src[k - x], c, _1 + _2*lenrek[k]) ;
               n += lenrek[k] ;
            }
            static_transform(c, *dst, _1 * kernel_sum/n) ;
         }

         // The bulk of the pixels in the middle of the image
         for (int x = 0; x < W - w + 1; ++x, ++src, ++dst)
         {
            pixel c(0) ;
            for (int k = 0; k < w; ++k)
               static_transform(c, src[k], c, _1 + _2*lenrek[k]) ;
            *dst = c ;
         }

         // Right-most w2 pixels
         for (int x = w - 1; x > w2; --x, ++src, ++dst)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = 0; k < x; ++k) {
               static_transform(c, src[k], c, _1 + _2*lenrek[k]) ;
               n += lenrek[k] ;
            }
            static_transform(c, *dst, _1 * kernel_sum/n) ;
         }
      }
   }
   return output_img ;
}

/**
   \brief  Convolve input image with supplied kernel along x-direction.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic xfilter
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type
xfilter(const image_type& img, const std::vector<float>& kernel)
{
   return xfilter<image_type, image_type>(img, kernel) ;
}

/**
   \brief  Convolve input image with supplied kernel along y-direction.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   This function applies the one-dimensional kernel (specified via an
   STL vector of floats) to the input image along the y-direction and
   returns a new image object containing the filteration result.

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename output_image_type, typename input_image_type>
input_image_type
yfilter(const input_image_type& img, const std::vector<float>& kernel)
{
   namespace gil = boost::gil ;

   const int W  = img.width() ;
   const int H  = img.height();
   const int h  = kernel.size() ;
   const int h2 = (h - 1)/2 ;

   // Flip the kernel (this makes the convolution a little more efficient)
   // DEVNOTE: lenrek = kernel spelled backwards, i.e., flipped kernel.
   std::vector<float> lenrek(h) ;
   float kernel_sum = 0 ;
   for (int i = 0; i < h; ++i) {
      kernel_sum += kernel[i] ;
      lenrek[h - 1 - i] = kernel[i] ;
   }

   typedef typename input_image_type::const_view_t in_view ;
   typedef typename in_view::x_iterator in_iterator ;
   typedef typename in_view::value_type pixel ;

   typedef typename output_image_type::view_t out_view ;
   typedef typename out_view::x_iterator out_iterator ;

   // Create a result image and setup "pointers" into the source and
   // destination pixel arrays.
   output_image_type output_img(W, H) ;
   in_view  in  = gil::const_view(img) ;
   out_view out = gil::view(output_img);

   using gil::static_transform ;
   using boost::lambda::_1 ;
   using boost::lambda::_2 ;

   // Do the convolution...
   if (H < h) // special case for small images (i.e., smaller than the kernel)
   {
      for (int y = 0; y < H; ++y)
      {
         in_iterator  src =  in.row_begin(y) ;
         out_iterator dst = out.row_begin(y) ;
         for (int x = 0; x < W; ++x, ++src, ++dst)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = 0; k < h; ++k)
               if (y + k - h2 >= 0 && y + k - h2 < H) {
                  static_transform(c, src[W * (k - h2)], c, _1 + _2*lenrek[k]) ;
                  n += lenrek[k] ;
               }
            static_transform(c, *dst, _1 * kernel_sum/n) ;
         }
      }
   }
   else // convolution of regular sized images (i.e., larger than the kernel)
   {
      // Since we're moving vertically along an image stored in row-major
      // order, consecutive elements will have their indices every W
      // spaces apart. Rather than compute these indices over and over
      // inside the convolution loops, we precompute them.
      std::vector<int> index(h) ;
      for (int i = 0; i < h; ++i)
         index[i] = W * i ;

      int row = 0 ;

      // Top-most h2 pixels
      for (int y = h2; y > 0; --y, ++row)
      {
         in_iterator  src =  in.row_begin(0) ;
         out_iterator dst = out.row_begin(row) ;
         for (int x = 0; x < W; ++x, ++src, ++dst)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = y; k < h; ++k) {
               static_transform(c, src[index[k - y]], c, _1 + _2*lenrek[k]) ;
               n += lenrek[k] ;
            }
            static_transform(c, *dst, _1 * kernel_sum/n) ;
         }
      }

      // The bulk of the pixels in the middle of the image
      for (int y = 0; y < H - h + 1; ++y, ++row)
      {
         in_iterator  src =  in.row_begin(row - h2) ;
         out_iterator dst = out.row_begin(row) ;
         for (int x = 0; x < W; ++x, ++src, ++dst)
         {
            pixel c(0) ;
            for (int k = 0; k < h; ++k)
               static_transform(c, src[index[k]], c, _1 + _2*lenrek[k]) ;
            *dst = c ;
         }
      }

      // The bottom-most h2 pixels
      for (int y = h - 1; y > h2; --y, ++row)
      {
         in_iterator  src =  in.row_begin(row - h2) ;
         out_iterator dst = out.row_begin(row) ;
         for (int x = 0; x < W; ++x, ++src, ++dst)
         {
            pixel c(0) ;
            float n(0) ;
            for (int k = 0; k < y; ++k) {
               static_transform(c, src[index[k]], c, _1 + _2*lenrek[k]) ;
               n += lenrek[k] ;
            }
            static_transform(c, *dst, _1 * kernel_sum/n) ;
         }
      }
   }
   return output_img ;
}

/**
   \brief  Convolve input image with supplied kernel along y-direction.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic yfilter
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type
yfilter(const image_type& img, const std::vector<float>& kernel)
{
   return yfilter<image_type, image_type>(img, kernel) ;
}

/**
   \brief  Convolve input image with supplied kernel.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   This function convolves the input image with the supplied kernel. The
   convolution operation is implemented as a separable filter wherein the
   input image is first convolved with the kernel along the x-direction.
   The result of this convolution is then convolved with the kernel along
   the y-direction.

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename out, typename in>
inline out convolve(const in& img, const std::vector<float>& kernel)
{
   return yfilter<out, in>(xfilter<out, in>(img, kernel), kernel) ;
}

/**
   \brief  Convolve input image with supplied kernel.
   \param  img Reference to the input image to be convolved.
   \param  kernel Reference to a one-dimensional kernel.
   \return The filtered image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic convolve
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type
convolve(const image_type& img, const std::vector<float>& kernel)
{
   return convolve<image_type, image_type>(img, kernel) ;
}

/**
   \brief  Blur image with Gaussian kernel.
   \param  img The input image to be blurred.
   \param  sigma The standard deviation for the Gaussian kernel.
   \param  kernel_width Half-size of the Gaussian kernel.
   \return Blurred image.

   This function applies a Gaussian blur to the input image and returns a
   new image object containing the blurred image.

   By default, the standard deviation for the Gaussian is 1.6 and the
   kernel half-size 4 (making for a total kernel size of 4*2 + 1 = 9).

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename out, typename in>
inline out blur(const in& img, float sigma = 1.6f, int kernel_width = 4)
{
   std::vector<float> gaussian_kernel(float, int) ; // function prototype
   return convolve<out, in>(img, gaussian_kernel(sigma, kernel_width)) ;
}

/**
   \brief  Blur image with Gaussian kernel.
   \param  img The input image to be blurred.
   \param  sigma The standard deviation for the Gaussian kernel.
   \param  kernel_width Half-size of the Gaussian kernel.
   \return Blurred image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic blur
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type
blur(const image_type& img, float sigma = 1.6f, int kernel_width = 4)
{
   return blur<image_type, image_type>(img, sigma, kernel_width) ;
}

//---------------------------- DECIMATION -------------------------------

/**
   \brief Reduce size of input image by given x and y factors.
   \param img Input image that is to be decimated.
   \param X Factor by which width should be scaled.
   \param Y Factor by which height should be scaled.
   \return New image object containing decimated image.

   This function scales the input image by the specified factors
   (default: halve both width and height).

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename output_image_type, typename input_image_type>
output_image_type
decimate(const input_image_type& img, int X = 2, int Y = 2)
{
   namespace gil = boost::gil ;

   // Don't go smaller than 1x1
   if (img.width() <= 1 && img.height() <= 1)
      return img ;

   if (img.width() == 1) // only need to thin out vertically
   {
      const int h = img.height()/Y;
      output_image_type output_img(1, h) ;

      typename  input_image_type::const_view_t src = gil::const_view(img) ;
      typename output_image_type::view_t       dst = gil::view(output_img);

      int s = 0, d = 0 ;
      for (int y = 0; y < h; ++y)
         dst[d++] = src[s++ + y * Y] ;
      return output_img ;
   }

   if (img.height() == 1) // only need to thin out horizontally
   {
      const int w = img.width()/X ;
      output_image_type output_img(w, 1) ;

      typename  input_image_type::const_view_t src = gil::const_view(img) ;
      typename output_image_type::view_t       dst = gil::view(output_img);

      int s = 0, d = 0 ;
      for (int x = 0; x < w; ++x)
         dst[d++] = src[s++ + x * X] ;
      return output_img ;
   }

   // Otherwise: decimate "regular" image (i.e., image whose dimensions
   // are greater than 1x1, 1xh or wx1).
   const int w = img.width ()/X ;
   const int h = img.height()/Y ;
   const int skip = img.width() % X + img.width() * (Y - 1) ;

   output_image_type output_img(w, h) ;
   typename  input_image_type::const_view_t src = gil::const_view(img) ;
   typename output_image_type::view_t       dst = gil::view(output_img);

   int s = 0, d = 0 ;
   for (int y = 0; y < h; ++y)
   {
      for (int x = 0; x < w; ++x)
      {
         dst[d++] = src[s] ;
         s += X ;
      }
      s += skip ;
   }
   return output_img ;
}

/**
   \brief Reduce size of input image by given x and y factors.
   \param img Input image that is to be decimated.
   \param X Factor by which width should be scaled.
   \param Y Factor by which height should be scaled.
   \return New image object containing decimated image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic decimate
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type decimate(const image_type& img, int X = 2, int Y = 2)
{
   return decimate<image_type, image_type>(img, X, Y) ;
}

//-------------------------- IMAGE DENOISING ----------------------------

/**
   \brief  Remove salt-n-pepper noise from input image.
   \param  img The input image that is to be denoised.
   \return New image object containing the denoised image.

   This function applies linear interpolation followed by a Gaussian blur
   followed by a decimation back to the original size in order to remove
   salt-n-pepper noise from the input image.

   As this is a template function, it can work with any image type
   supported by Boost.GIL. Note that we specify the output image type
   first and then the input image type in the template's type list. This
   allows clients to get away with specifying only the output image type
   and relying on type deduction for the input image type.
*/
template<typename out, typename in>
inline out denoise(const in& img)
{
   return decimate<out, in>(blur<out, in>(interpolate<out, in>(img))) ;
}

/**
   \brief  Remove salt-n-pepper noise from input image.
   \param  img The input image that is to be denoised.
   \return New image object containing the denoised image.

   Often, the input and output image types are the same. This is a
   convenience function that allows clients to call the generic denoise
   function with minimum fuss (using parameter type deduction) so that
   the call site does not have to be cluttered with unnecessary template
   verbosity.
*/
template<typename image_type>
inline image_type denoise(const image_type& img)
{
   return denoise<image_type, image_type>(img) ;
}

//-------------------------- IMAGE GRADIENTS ----------------------------

/// This type specifies a pixel that is a 2D vector that can be used to
/// represent a gradient image.
/// \return
typedef boost::gil::pixel<float,
                          boost::gil::layout<boost::gil::devicen_t<2> > >
        gradient_pixel_t ;

/// This type specifies an image of 2D vectors that can be used to
/// represent a gradient image.
/// \return
typedef boost::gil::image<gradient_pixel_t, false> gradient_image_t ;

/**
   \brief  Calculate the gradient of an image.
   \param  img The input image whose gradient we want to calculate.
   \return An image of 2D vectors containing the gradient of the input image.

   This function computes the gradient vector field of the input image
   that shows, at each pixel location, the magnitude and direction of
   maximal brightness. It returns this vector field using a Boost.GIL
   image whose pixels are 2D vectors (in polar form) indicating the
   gradient.
*/
template<typename image_type>
gradient_image_t gradient(const image_type& img)
{
   namespace gil = boost::gil ;

   typedef typename image_type::const_view_t in_view ;
   typedef typename in_view::xy_locator in_loc ;
   typedef typename in_loc::cached_location_t in_loc_cached ;
   typedef typename gradient_image_t::view_t out_view;
   typedef typename out_view::xy_locator out_loc ;

   const int w = img.width() ;
   const int h = img.height();

   in_view in  = gil::const_view(img) ;
   in_loc  src = in.xy_at(0,1) ;
   in_loc_cached src_e = src.cache_location(+1, 0) ;
   in_loc_cached src_w = src.cache_location(-1, 0) ;
   in_loc_cached src_n = src.cache_location( 0,-1) ;
   in_loc_cached src_s = src.cache_location( 0,+1) ;

   gradient_image_t g(w, h) ;
   out_view out = gil::view(g) ;
   out_loc  dst = out.xy_at(0,0) ;

   // First row of gradient vector field (GVF) is all zeros
   gradient_pixel_t zero(0) ;
   std::fill_n(dst.x(), w, zero) ;

   // Compute the GVF for the inner rows
   dst = out.xy_at(0,1) ;
   for (int y = 1; y < h - 1; ++y)
   {
      // Left-most pixel is zero
      *dst = zero ;
      ++src.x() ;
      ++dst.x() ;

      // Compute GVF for inner columns
      for (int x = 1; x < w - 1; ++x, ++src.x(), ++dst.x())
      {
         float dx = src[src_e] - src[src_w] ;
         float dy = src[src_s] - src[src_n] ;
         *dst = gradient_pixel_t(sqrtf(dx*dx + dy*dy), atan2f(dy, dx)) ;
      }

      // Right-most pixel is zero
      *dst = zero ;
      src += gil::point2<std::ptrdiff_t>(-w + 1, 1) ;
      dst += gil::point2<std::ptrdiff_t>(-w + 1, 1) ;
   }

   // The last row of the GVF is all zeros
   std::fill_n(dst.x(), w, zero) ;

   return g ;
}

//-----------------------------------------------------------------------

} // end of namespace encapsulating this file's definitions

#endif

/* So things look consistent in everyone's emacs... */
/* Local Variables: */
/* indent-tabs-mode: nil */
/* End: */
